%matplotlib inline

import numpy as np
import pandas as pd
from mpl_toolkits import mplot3d
import matplotlib.pyplot as plt
import scipy.stats as stats
from scipy.stats import pearsonr
import scipy.signal as signal
import seaborn as sns
import scipy.fft

G = 9.80665
SOURCE = 'esp'
# DATASET 'L3_StnVinohrady_Riazanska_Skoda30T.csv'
# DATASET = 'L3_Pionierska_RacianskeMyto_Skoda30T.csv'
# DATASET = 'L4_7954_ZaluhyKrizovatka_KutikyObratisko.csv'
# DATASET = 'L9_Postova_KralovskeUdolie_Skoda29T.csv'
# DATASET = 'L9_Lanfranconi_Riviera_Skoda29T.csv'
# DATASET = 'L20_3014_Zaluhy-Drobneho.csv'
# DATASET = 'L35_1915_Most_Kutiky_neutral.csv'
# DATASET = 'L35_1915_Borska_Zaluhy.csv'
# DATASET = 'L83_4940_PriKrizi_Alexyho.csv'
# DATASET = 'L83_4940_Alexyho_Svantnerova.csv'
 
SOURCE = 'r-das'
# DATASET = 'empty_motor_952hz.txt'
# DATASET = 'motor_with_small_cd_952hz.txt'
# DATASET = 'motor_with_small_cd_unbalanced_952hz.txt'
# DATASET = 'servo952hz.txt'
# DATASET = 'smd_vibrator_2.35V_2g_952hz_held.txt'
# DATASET = 'smd_vibrator_2.35V_4g_952hz_held.txt'
# DATASET = 'smd_vibrator_2.35V_4g_952hz.txt'
# DATASET = 'smd_vibrator_2.45V_952hz.txt'
# DATASET = 'st_meranie_s_hybatkom_2hz_952hz.txt'
# DATASET = 'st_meranie_s_hybatkom_2hz_na_lukasa_952hz.txt'
# DATASET = 'st_meranie_s_hybatkom_4hz_952hz.txt'
DATASET = 'st_meranie_s_hybatkom_4hz_na_lukasa_952hz.txt'

# Výstupom je tabuľka s troma atribútmi osí akcelerácie [m/s^2]: x, y, z a časovým indexom t [s]
if SOURCE == 'r-das': 
    FS_HZ = 952
    PATH = 'datasets/r-das/'
    
    measurement = pd.read_csv(PATH + DATASET, sep='\t', skiprows=2)
    measurement = measurement[['A_X [mg] ', ' A_Y [mg] ', ' A_Z [mg] ']].rename(columns={
        'A_X [mg] ': 'x', ' A_Y [mg] ': 'y', ' A_Z [mg] ': 'z'
    })
    measurement['t'] = measurement.index * (1 / FS_HZ)
    measurement.set_index('t', inplace=True)
    vibrations = measurement.transform(lambda v: (v / 1000) * G)   # Premena z mg na m/s^2
    
elif SOURCE == 'esp':
    FS_HZ = 500
    MG_LSB = 0.061
    PATH = 'datasets/esp/'
    
    measurement = pd.read_csv(PATH + DATASET, sep=' ').rename(columns={
        'A_X[LSB]': 'x', 'A_Y[LSB]': 'y', 'A_Z[LSB]': 'z'
    })
    measurement['t'] = measurement.index * (1 / FS_HZ)
    measurement.set_index('t', inplace=True)
    vibrations = measurement.transform(lambda v: ((v * MG_LSB) / 1000) * G)   # Premena z LSB na m/s^2
        
    # Odstrániť chybné hodnoty mimo rozsah: 
    # vibrations[vibrations['y'].abs() > 2*G]
    # measurement[measurement.isnull().any(axis=1)]
    # pd.to_numeric(measurement['x'])

else:
    raise NotImplementedError('Neplatný názov datasetu')

vibrations.head()

	x	y	z
t
0.000000	0.039227	3.285228	5.011198
0.001050	-2.147656	9.414384	7.933580
0.002101	-6.364516	19.064128	14.739395
0.003151	-1.706357	7.354988	11.346294
0.004202	-0.794339	-10.954028	6.796008

vibrations.info()

<class 'pandas.core.frame.DataFrame'>
Float64Index: 14812 entries, 0.0 to 15.557773109243696
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       14812 non-null  float64
 1   y       14812 non-null  float64
 2   z       14812 non-null  float64
dtypes: float64(3)
memory usage: 462.9 KB

vibrations.describe()

	x	y	z
count	14812.000000	14812.000000	14812.000000
mean	-0.285808	0.441458	9.973866
std	4.271749	7.113732	3.737934
min	-19.593687	-19.593687	-4.501252
25%	-2.520309	-4.265893	7.531507
50%	-0.215746	0.568786	9.885103
75%	1.951523	5.060231	12.346572
max	16.926278	19.593687	19.593687

pd.DataFrame([
    ['mean(A):', *[vibrations[c].mean() for c in vibrations.columns.values]],
    ['stdev(A):', *[vibrations[c].std() for c in vibrations.columns.values]],
    ['skew(A):', *[stats.skew(vibrations[c]) for c in vibrations.columns.values]],
    ['kurtosis(A):', *[stats.kurtosis(vibrations[c]) for c in vibrations.columns.values]],
    ['MAD(A):', *[stats.median_abs_deviation(vibrations[c]) for c in vibrations.columns.values]]
], columns=['Moment', 'x', 'y', 'z']).set_index('Moment')

	x	y	z
Moment
mean(A):	-0.285808	0.441458	9.973866
stdev(A):	4.271749	7.113732	3.737934
skew(A):	-0.033331	-0.002547	0.009820
kurtosis(A):	1.031522	-0.183686	0.166633
MAD(A):	2.235916	4.697385	2.402629

sns.set(rc={'figure.figsize': (10, 6)})
sns.boxplot(x='variable', y='value', data=pd.melt(vibrations))

<AxesSubplot:xlabel='variable', ylabel='value'>

sns.set(rc={'figure.figsize': (20, 8)})
df = vibrations.reset_index().melt('t', var_name='cols',  value_name='vals')
sns.lineplot(x='t', y='vals', hue='cols', data=df)

<AxesSubplot:xlabel='t', ylabel='vals'>

Amplitúda akcelerácie

sns.set(rc={'figure.figsize': (20, 12)})
mag = np.sqrt(vibrations.x**2 + vibrations.y**2 + vibrations.z**2)

fig, axes = plt.subplots(3, 1)

axes[0].set_title('Amplitúda')
axes[1].set_title('Amplitúda korigovaná')
axes[2].set_title('Magnitúda')

sns.lineplot(x='t', y=mag, data=vibrations, ax=axes[0])
sns.lineplot(x='t', y=mag - mag.mean(), data=vibrations, ax=axes[1])
sns.lineplot(x='t', y=np.abs(mag - mag.mean()), data=vibrations, ax=axes[2])

<AxesSubplot:title={'center':'Magnitúda'}, xlabel='t'>

Časový priebeh

plt.style.use('default')
plt.rcParams['figure.figsize'] = (20, 15)
vibrations.plot(subplots=True, grid=True)

array([<AxesSubplot:xlabel='t'>, <AxesSubplot:xlabel='t'>,
       <AxesSubplot:xlabel='t'>], dtype=object)

Vyhladzovanie filtrom cez konvolučné jadro

N = 32
kernel = np.ones(N) / N

smooth = pd.DataFrame({
    't': vibrations.index.values[kernel.size - 1:], 
    'x': np.convolve(vibrations['x'].values, kernel, mode='valid'),
    'y': np.convolve(vibrations['y'].values, kernel, mode='valid'),
    'z': np.convolve(vibrations['z'].values, kernel, mode='valid')
}).set_index('t')

plt.style.use('seaborn')
plt.rcParams['figure.figsize'] = (20, 12)
smooth.plot(subplots=True)

array([<AxesSubplot:xlabel='t'>, <AxesSubplot:xlabel='t'>,
       <AxesSubplot:xlabel='t'>], dtype=object)

Orezanie hodnôt s nízkou magnitúdou (efekt vzorkovania iba pri zdroji prerušenia)

THETA = 0.5

trim = pd.DataFrame({
    't': vibrations.index,
    'x': vibrations['x'],
    'x_event': vibrations['x'].where(vibrations['x'].abs() > THETA, 0)
}).set_index('t')

plt.style.use('seaborn')
plt.rcParams['figure.figsize'] = (20, 10)
trim.plot(subplots=True)

print('Redukcia dát o %.4f %%' % ((1 - (trim['x_event'][trim['x_event'] != 0].size / trim['x'].size)) * 100))
plt.show()

Redukcia dát o 13.6984 %

x = np.linspace(0, 12)
y = pd.Series(np.sin(x))
y_cut = y.where(y.abs() > THETA, 0)

plt.rcParams['figure.figsize'] = (20, 5)
plt.plot(x, y)
plt.plot(x, y_cut)

[<matplotlib.lines.Line2D at 0x7f018d193f40>]

n = trim['x'].size
x = np.fft.fftfreq(n, 1 / FS_HZ)[:n//2]

y = np.absolute(np.fft.fft(trim['x']))[:n//2]
y = 20 * np.log10(y / y.max())

yf = np.absolute(np.fft.fft(trim['x_event']))[:n//2] 
yf = 20 * np.log10(yf / yf.max())

plt.style.use('seaborn')
plt.rcParams['figure.figsize'] = (20, 5)
plt.xlim(0, FS_HZ // 2)
plt.plot(x, y, label='x', alpha=0.9)
plt.plot(x, yf, label='x_event', alpha=0.9)

[<matplotlib.lines.Line2D at 0x7f018d112d60>]

Štatistické momenty podľa okien

plt.style.use('default')
plt.rcParams['figure.figsize'] = (20, 15)

W = 512                   # Veľkosť okna
OVERLAP = 0               # Veľkosť prekrytia okien <0; W)
STEP = W - OVERLAP        # Posun okna na základe prekrytia
AXIS = 'x'

fig, axes = plt.subplots(3, 1)
x = vibrations[AXIS].values
t = [vibrations.index[i] for i in range(0, len(x) - W, STEP)]

for i, col in enumerate(vibrations.columns.values, start=0):
    x = vibrations[col].values
    axes[i].plot(t, [np.mean(x[i:i+W]) for i in range(0, len(x) - W, STEP)], label='mean', marker='o')
    axes[i].plot(t, [np.median(x[i:i+W]) for i in range(0, len(x) - W, STEP)], label='median', marker='o')
    axes[i].plot(t, [np.std(x[i:i+W]) for i in range(0, len(x) - W, STEP)], label='std', marker='o')
    axes[i].plot(t, [stats.skew(x[i:i+W]) for i in range(0, len(x) - W, STEP)], label='skew', marker='o')
    axes[i].plot(t, [stats.kurtosis(x[i:i+W]) for i in range(0, len(x) - W, STEP)], label='kurtosis', marker='o')
    axes[i].plot(t, [stats.median_abs_deviation(x[i:i+W]) for i in range(0, len(x) - W, STEP)], label='MAD', marker='o')
    axes[i].legend()
    axes[i].set_xlabel('Čas [s]')
    axes[i].set_ylabel('Zrýchlenie [m/s^2]')
    axes[i].grid()

Korelácie podľa okien

Z dokumentácie (https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.pearsonr.html):

• Correlations of -1 or +1 imply an exact linear relationship. Positive correlations imply that as x increases, so does y. Negative correlations imply that as x increases, y decreases.

• The p-value roughly indicates the probability of an uncorrelated system producing datasets that have a Pearson correlation at least as extreme as the one computed from these datasets. (Čiže p < alpha znamená že sú korelované)

plt.style.use('default')
plt.rcParams['figure.figsize'] = (20, 10)

t = [vibrations.index[i] for i in range(0, len(x) - W, STEP)]
x = vibrations['x'].values
y = vibrations['y'].values
z = vibrations['z'].values
coor_xy = [pearsonr(x[i:i+W], y[i:i+W]) for i in range(0, len(x) - W, STEP)]
coor_xz = [pearsonr(x[i:i+W], z[i:i+W]) for i in range(0, len(x) - W, STEP)]
coor_yz = [pearsonr(y[i:i+W], z[i:i+W]) for i in range(0, len(x) - W, STEP)]

fig, axes = plt.subplots(2, 1)
for i in range(2):
    axes[i].plot(t, [n[i] for n in coor_xy], label='XY (Yaw)', marker='o')
    axes[i].plot(t, [n[i] for n in coor_xz], label='XZ (Pitch)', marker='o')
    axes[i].plot(t, [n[i] for n in coor_yz], label='YZ (Roll)', marker='o')
    axes[i].grid(True)
    axes[i].legend()
    axes[i].set_xlabel('Čas [s]')

axes[0].set_ylabel('Osová korelácia')
axes[1].set_ylabel('p-hodnota korelácie')   # Nižšie ako <0.05 je korelované
plt.show()

Histogram rozdelení zložiek akcelerácie

plt.style.use('default')
plt.rcParams['figure.figsize'] = (15, 5)
bins = 200
plt.hist(vibrations['x'] - vibrations['x'].mean(), bins, alpha=0.5, label='ax', edgecolor='k')
plt.hist(vibrations['y'] - vibrations['y'].mean(), bins, alpha=0.5, label='ay', edgecolor='k')
plt.hist(vibrations['z'] - vibrations['z'].mean(), bins, alpha=0.5, label='az', edgecolor='k')
plt.legend()
plt.show()

Detekcia anomálií Z-skóre algoritmom

plt.rcParams['figure.figsize'] = (20, 6)
K = 3
AXIS = 'y'
r = vibrations[abs(vibrations[AXIS] - vibrations[AXIS].mean()) >= K * vibrations[AXIS].std()][AXIS]

ax = vibrations[AXIS].plot()
ax.plot(r.index.values, r.values, 'o', color='r')
ax.set_xlabel('Čas [s]')
p = ax.set_ylabel('Zrýchlenie [m/s^2]')

Spektrogram

W = 256
%matplotlib inline
plt.style.use('default')
plt.rcParams['figure.figsize'] = (20, 15)
fig, axes = plt.subplots(3, 1)

for i, col in enumerate(vibrations.columns.values, start=0):
    f, t, Sxx = signal.spectrogram(
        vibrations[col], FS_HZ, nperseg=W, mode='magnitude', window='tukey', noverlap=W//8, scaling='spectrum'
    )
    axes[i].pcolormesh(t, f, Sxx, shading='gouraud')
    axes[i].set_ylabel('Frekvencia [Hz] ' + col)
    # axes[i].set_ylim(0, 50)
q = plt.xlabel('Čas [s]')

Logaritmický spektrogram s nastaviteľným vzájomným prekryvom okien

plt.style.use('default')
plt.rcParams['figure.figsize'] = (20, 5)

W = 512                 # Veľkosť okna
OVERLAP = 0            # Veľkosť prekrytia okien <0; W)
STEP = W - OVERLAP       # Posun okna na základe prekrytia
window = signal.tukey(W) # Window function. Pre boxcar: np.ones(W)
FS = 1/FS_HZ             # Vzorkovacia frekvencia

t = vibrations.index
x = vibrations['x'].values

t = np.array([t[i] for i in range(0,  len(x) - W, STEP)])
f = np.fft.fftfreq(W, FS)[:W//2]
f_mag = np.asarray([
    np.log(np.absolute(np.fft.fft(x[i:i+W] * window)[:W//2]))
    for i in range(0, len(x) - W, STEP)
]).T

plt.pcolormesh(t, f, f_mag, shading='gouraud')
plt.ylabel('Frekvencia [Hz]')
plt.xlabel('Čas [s]')
plt.show()

Porovnanie DCT

plt.style.use('default')
plt.rcParams['figure.figsize'] = (20, 30)

W = 512                   # Veľkosť okna
OVERLAP = 0               # Veľkosť prekrytia okien <0; W)
STEP = W - OVERLAP        # Posun okna na základe prekrytia
window = signal.tukey(W)  # Window function. Pre boxcar: np.ones(W)
FS = 1 / FS_HZ            # Vzorkovacia frekvencia je 200 Hz

t = vibrations.index
x = vibrations['x'].values

t = np.array([t[i] for i in range(0,  len(x) - W, STEP)])
f = np.linspace(0, 1/FS/2, W)

# Linear
f_mag = []
for kind in range(1, 5):
    f_mag.append(np.asarray([
        np.abs(scipy.fft.dct(x[i:i+W], type=kind))
        for i in range(0, len(x) - W, STEP)
    ]).T)

# Logarithmic
f_mag_log = []
for kind in range(1, 5):
    f_mag_log.append(np.asarray([
        np.log(np.abs(scipy.fft.dct(x[i:i+W], type=kind)))
        for i in range(0, len(x) - W, STEP)
    ]).T)

    
fig, axes = plt.subplots(8, 1)
for i in range(len(f_mag)):
    axes[i].set_title(label='DCT type ' +  str(i + 1))
    axes[i].pcolormesh(t, f, f_mag[i], shading='gouraud')
    axes[i].set_ylabel('Frekvencia [Hz]')

for i in range(len(f_mag_log)):
    j = i + 4
    axes[j].set_title(label='DCT type ' +  str(i + 1))
    axes[j].pcolormesh(t, f, f_mag_log[i], shading='gouraud')
    axes[j].set_ylabel('Frekvencia [Hz]')

plt.xlabel('Čas [s]')
plt.show()

plt.rcParams['figure.figsize'] = (20, 12)
fig, axes = plt.subplots(3, 1)

for i, col in enumerate(vibrations.columns.values, start=0):
    f, t, Sxx = signal.spectrogram(vibrations[col], FS_HZ, nperseg=W, mode='magnitude', window='tukey', noverlap=W//8, scaling='spectrum')
    peaks = np.argmax(Sxx, axis=0)
    peaks = np.take(f, peaks)
    axes[i].scatter(t, peaks, color='m')
    axes[i].set_ylabel('Frekvencia [Hz]')
    axes[i].grid(True)

q = plt.xlabel('Čas [s]')

Frekvenčné zložky v okne (s oknovou funkciou) veľkosti N vzoriek s posunom T sekúnd v lineárnej a logaritmickej mierke

plt.style.use('default')
plt.rcParams['figure.figsize'] = (20, 12)

T = 8
N = 256
AXIS = 'z'
window = np.hamming(N)

plt.rcParams['figure.figsize'] = (20, 10)           
y = vibrations[AXIS][T:].values[0:N] * window

spectrum = np.fft.rfft(y)                               
magnitude = 2 * np.abs(spectrum) / np.sum(window)
# dbfs = 10 * np.log10(magnitude)
dbfs = 20 * np.log10(magnitude / magnitude.max())
frequencies = np.arange((N / 2) + 1) / (float(N) / FS_HZ)

fig, axes = plt.subplots(2, 1)

for i, col in enumerate((magnitude, dbfs), start=0):
    axes[i].plot(frequencies, col)
    axes[i].set_xlabel('Frekvencia [Hz]')
    axes[i].set_ylabel('Amplitúda')
    axes[i].set_ylim(col.min(), col.max())
    #axes[i].set_xscale('log')
    axes[i].grid(True)

Nájdenie špičiek

spikes, _ = signal.find_peaks(dbfs, height=-20, distance=5)

plt.rcParams['figure.figsize'] = (20, 5)
plt.plot(frequencies, dbfs, alpha=0.5)
plt.grid()

if spikes.size > 0:
    plt.plot(frequencies[spikes], dbfs[spikes], 'o', color='r')

pd.DataFrame({
    'Frekvencia [Hz]': frequencies[spikes],
    'Amplitútda [dB]': dbfs[spikes],
}).set_index('Frekvencia [Hz]').sort_values(ascending=False, by='Amplitútda [dB]').head(10)

	Amplitútda [dB]
Frekvencia [Hz]
171.0625	-19.72462

Numerická integrácia

%matplotlib inline
plt.style.use('seaborn')
plt.rcParams['figure.figsize'] = (20, 10)

motion = vibrations.reset_index().rename(columns={'x': 'ax', 'y': 'ay', 'z': 'az'})
motion['dt'] = motion['t'] - motion['t'].shift(-1)
motion.set_index('t', inplace=True)

# Obdĺžníkové pravidlo
motion['vx'] = np.cumsum(motion['ax'] * motion['dt'])
motion['vy'] = np.cumsum(motion['ay'] * motion['dt'])
motion['vz'] = np.cumsum(motion['az'] * motion['dt'])

# Lichobežníkové pravidlo
motion['vx_trapez'] = np.cumsum(((motion['ax'] + motion['ax'].shift(-1)) / 2) * motion['dt'])
motion['vy_trapez'] = np.cumsum(((motion['ay'] + motion['ay'].shift(-1)) / 2) * motion['dt'])
motion['vz_trapez'] = np.cumsum(((motion['ay'] + motion['az'].shift(-1)) / 2) * motion['dt'])

fig, ax = plt.subplots(3, 1)
motion['vx'].plot(ax=ax[0])
motion['vx_trapez'].plot(ax=ax[0])
ax[0].legend()

motion['vy'].plot(ax=ax[1])
motion['vy_trapez'].plot(ax=ax[1])

motion['vz'].plot(ax=ax[2])
motion['vz_trapez'].plot(ax=ax[2])

<AxesSubplot:xlabel='t'>

3D graf párov zložiek akcelerácie v čase

%matplotlib inline
plt.style.use('default')
plt.rcParams['figure.figsize'] = (20, 20)

fig = plt.figure()
ax = fig.add_subplot(2, 2, 1, projection='3d')
ax.scatter3D(vibrations.x, vibrations.y, vibrations.z, c=vibrations.z)

ax = fig.add_subplot(2, 2, 2, projection='3d')
ax.scatter3D(vibrations.index, vibrations.x, vibrations.y, c=vibrations.index)

ax = fig.add_subplot(2, 2, 3, projection='3d')
ax.scatter3D(vibrations.index, vibrations.x, vibrations.z, c=vibrations.index)

ax = fig.add_subplot(2, 2, 4, projection='3d')
ax.scatter3D(vibrations.index, vibrations.y, vibrations.z, c=vibrations.index)

plt.show()

Interaktívne frekvenčné spektrum

%matplotlib widget
import ipywidgets as widgets
from ipywidgets import Layout

plt.rcParams['figure.figsize'] = (10, 5) 


def spectrum(N=1024, T=0, WF='hamming', AXIS='x', TRANSFORM='dct-ii'):
    if WF == 'boxcar':
        window = np.ones(N)
    elif WF == 'hamming':
        window = np.hamming(N)
    elif WF == 'hanning':
        window = np.hanning(N)
    elif WF == 'blackman':
        window = np.blackman(N)
    elif WF == 'bartlett':
        window = np.bartlett(N)
        
    col = vibrations[AXIS]
    y = col[T:].values[0:N] * window
    
    if TRANSFORM == 'fft':
        spectrum = np.fft.fft(y)[:N//2+1]
    elif TRANSFORM == 'rfft':
        spectrum = np.fft.rfft(y)   
    elif TRANSFORM == 'dct-ii':
        spectrum = scipy.fft.dct(y, type=2)[:N//2+1]
    elif TRANSFORM == 'dct-iv':
        spectrum = scipy.fft.dct(y, type=4)[:N//2+1]
                                 
    magnitude = 2 * np.abs(spectrum) / np.sum(window)
    dbfs = 20 * np.log10(magnitude / magnitude.max())
    frequencies = np.arange((N / 2) + 1) / (float(N) / FS_HZ)
    
    return frequencies, dbfs


frequencies, dbfs = spectrum()
fig, axes = plt.subplots()

line, = axes.plot(frequencies, dbfs, marker='o', markersize=4)
axes.set_xlabel('Frekvencia [Hz]')
axes.set_ylabel('Amplitúda [dB]')
axes.set_ylim(-60, 5)
axes.grid(True)


window_slider = widgets.FloatLogSlider(
    value=512, base=2, min=3, max=12, step=1, 
    description='Veľkosť okna', 
    layout=Layout(width='50%'),
    continuous_update=True
)
t_slider = widgets.FloatSlider(
    value=1, min=0, max=t.size, step=0.1, 
    description='Čas [s]', 
    layout=Layout(width='50%'),
    continuous_update=True
)
transform_func = widgets.Dropdown(
    options=['fft', 'rfft', 'dct-ii', 'dct-iv'],
    value='rfft',
    description='Trasnformácia:',
    disabled=False,
)
win_func = widgets.Dropdown(
    options=['boxcar', 'hanning', 'hamming', 'blackman', 'bartlett'],
    value='hamming',
    description='Oknová funkcia:',
    disabled=False,
)
axis_choose = widgets.Dropdown(
    options=['x', 'y', 'z'],
    value='x',
    description='Os:',
    disabled=False,
)

def update(change):
    frequencies, dbfs = spectrum(
        N=int(window_slider.value), 
        T=int(t_slider.value), 
        WF=win_func.value,
        AXIS=axis_choose.value,
        TRANSFORM=transform_func.value
    )
    line.set_xdata(frequencies)
    line.set_ydata(dbfs)
    fig.canvas.draw()

    
window_slider.observe(update)
t_slider.observe(update)
transform_func.observe(update)
win_func.observe(update)
axis_choose.observe(update)

display(widgets.VBox(children=[transform_func, axis_choose, win_func, window_slider, t_slider]))

%matplotlib widget
import ipywidgets as widgets
from ipywidgets import Layout

plt.style.use('default')
plt.rcParams['figure.figsize'] = (9, 4) 


def spectrum(N=1024, OVERLAP=0, WF='hamming', AXIS='x'):
    if WF == 'boxcar':
        window = np.ones(N)
    elif WF == 'hamming':
        window = np.hamming(N)
    elif WF == 'hanning':
        window = np.hanning(N)
    elif WF == 'blackman':
        window = np.blackman(N)
    elif WF == 'bartlett':
        window = np.bartlett(N)
    
    STEP = N - OVERLAP
    t = vibrations.index
    x = vibrations[AXIS].values
                                 
    t = np.array([t[i] for i in range(0,  len(x) - N, STEP)])
    f = np.fft.fftfreq(N, 1 / FS_HZ)[:N//2]
    f_mag = np.asarray([
        np.log(np.absolute(np.fft.fft(x[i:i+N] * window)[:N//2]))
        for i in range(0, len(x) - N, STEP)
    ]).T

    return t, f, f_mag


t, f, f_mag = spectrum()
fig, axes = plt.subplots()

axes.pcolormesh(t, f, f_mag, shading='gouraud')
axes.set_ylabel('Frekvencia [Hz]')
axes.set_xlabel('Čas [s]')
axes.grid(True)


window_slider = widgets.FloatLogSlider(
    value=512, base=2, min=3, max=12, step=1, 
    description='Veľkosť okna', 
    layout=Layout(width='50%'),
    continuous_update=False
)
win_func = widgets.Dropdown(
    options=['boxcar', 'hanning', 'hamming', 'blackman', 'bartlett'],
    value='hamming',
    description='Oknová funkcia:',
    disabled=False,
)
axis_choose = widgets.Dropdown(
    options=['x', 'y', 'z'],
    value='x',
    description='Os:',
    disabled=False,
)

def update(change):
    t, f, f_mag = spectrum(
        N=int(window_slider.value),
        WF=win_func.value,
        AXIS=axis_choose.value,
    )
    axes.pcolormesh(t, f, f_mag, shading='gouraud')
    fig.canvas.draw()

    
window_slider.observe(update)
transform_func.observe(update)
win_func.observe(update)
axis_choose.observe(update)
t_slider.observe(update)

display(widgets.VBox(children=[axis_choose, win_func, window_slider]))